Analyze the ABC output

using CSV,DataFrames,Plots,Statistics,Distributions,Measures,StatsPlots,ColorSchemes
gr();

include("auxilliary.jl"); include("parameters.jl"); include("flow.jl"); include("abcmc.jl");
ENV["COLUMNS"]=200;

Inputs

# Indices of the ABC files that store the ABC parameter and trajectory files
frt = 1; lst = 10;

# threshold for ABC RMSE error cutoff, given as a top quantile (eg 0.05 is top 5%)
qℓ = 0.1;

Plot formatting

# Plot formatting
#  Finding that to get right output dimensions and font sizes for gr have to manually
#  compute scaling factors based on an observed print output
λplot = 100*7.5/10.42; λfont = 12/18; #scaling factors

# Final output two panel figure dimensions
figwdth = 7.45; fighght = 3.25; # inches

# Final output figure fontsizes
figtickfontsize = 12; figguidefontsize = 14; figtitlefontsize = 16;

lnopt = [:solid, :dash, :dot, :dashdot, :dashdotdot];
mkopt = [:none, :auto, :circle, :rect, :star5, :diamond, :hexagon, :cross, 
         :xcross, :utriangle, :dtriangle, :rtriangle, :ltriangle, :pentagon, 
         :heptagon, :octagon, :star4, :star6, :star7, :star8, :vline, :hline, 
         :+, :x];

Extract samples

df = CSV.read("ABCsmp$frt.csv",DataFrame,header=false);
for id=frt+1:lst
    dftemp = CSV.read("ABCsmp$id.csv",DataFrame,header=false);
    df = hcat(df,dftemp[:,2:end],makeunique=true);
end

# Dataframe of trajectories
dftrj = CSV.read("ABCtrj$frt.csv",DataFrame,header=false);
for id=frt+1:lst
    dftemp = CSV.read("ABCtrj$id.csv",DataFrame,header=false);
    dftrj = hcat(dftrj,dftemp[:,2:end],makeunique=true);
end

ncols = ncol(df)-1; nℓcols = qℓ*ncols;
println("The total number of abc sampling before conditioning: $ncols")
println("The total number of abc sampling after conditioning with top $qℓ quantile: $nℓcols")

The total number of abc sampling before conditioning: 5000
The total number of abc sampling after conditioning with top 0.1 quantile: 500.0

prmrg,prmvary=abcdata();prmvary[:ℓerr]=true;prmvary[:βθ]=true;
println("Varied parameters:")
for key in keys(prmvary)
    if prmvary[key]
        println(key)
    end
end

Varied parameters:
ρ
ℓerr
βθ
αeff
βα
γα
γθ

# Create dictionary of row position of key parameters;
mykeys = [key for key in keys(prmvary)]; mykeys=vcat(mykeys,[:ℓerr,:βθ])
pnt=Dict{Symbol,Int64}()
for key in mykeys
    pos = 1;
    while df[pos,1]!=string(key)
        pos+=1;
    end
    pnt[key]=pos;
end;

# Truncate down to trajectories that didn't fail to converge at their numerical tolerances
function Base.isnan(s::Union{Symbol,AbstractString})
    return false
end

dfnan = .!(isnan.(df[pnt[:ℓerr]:pnt[:ℓerr],:])); flagnan = [dfnan[1,i] for i=1:ncol(dfnan)];
df = df[:,flagnan]; dftrj = dftrj[:,flagnan];
print("The number of sims that failed to converge is "); println(sum(.!(flagnan)))

The number of sims that failed to converge is 1

# Extract marginal values conditioned on threshold
marginals = Dict{Symbol,Vector{Float64}}();
for key in keys(prmvary)
    if prmvary[key]
        marginals[key] = [v for v in df[pnt[key],2:end]];
    end
end
flag = marginals[:ℓerr] .<= quantile(marginals[:ℓerr],qℓ);
for key in keys(marginals)
    marginals[key] = marginals[key][flag]
end;

# Extract trajectories conditioned on threshold
dftrjabc = dftrj[:,2:end];
dftrjabc = dftrjabc[:,flag];
dftrjabc = hcat(DataFrame("day"=>convert(Vector,0:(nrow(dftrjabc)-1))),dftrjabc);
last(dftrjabc,7)

7 rows × 501 columns (omitted printing of 481 columns)

	day	Column3	Column14	Column17	Column22	Column30	Column34	Column39	Column57	Column61	Column63	Column66	Column72	Column85	Column96	Column105	Column106	Column107	Column117	Column130
	Int64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	55	8287.16	6054.01	8596.8	8241.08	8778.35	9934.67	9871.39	8675.37	8851.32	9589.75	9502.8	8200.3	7159.61	9212.67	8619.24	7965.41	9022.96	9555.09	8301.21
2	56	8276.62	6001.05	8493.72	8212.31	8766.59	10019.4	9918.44	8669.85	8775.63	9636.07	9557.39	8179.52	6934.02	9143.5	8614.97	7962.06	9077.26	9594.97	8123.06
3	57	8263.3	5950.36	8382.95	8171.35	8738.56	10095.1	9935.27	8655.34	8703.01	9655.54	9609.01	8149.9	6714.52	9062.82	8606.68	7956.75	9115.63	9612.36	7958.81
4	58	8252.61	5901.46	8253.94	8123.48	8700.8	10165.9	9917.77	8635.0	8617.94	9647.61	9661.05	8112.4	6486.68	8960.48	8597.13	7951.05	9125.09	9601.81	7800.75
5	59	8248.46	5853.46	8099.61	8074.53	8660.98	10236.0	9867.93	8612.29	8507.5	9616.81	9715.97	8068.97	6240.1	8829.81	8588.49	7946.12	9098.05	9563.43	7640.3
6	60	8252.28	5805.3	7917.44	8029.44	8625.77	10308.9	9793.02	8590.14	8363.58	9573.7	9774.57	8021.95	5969.88	8668.71	8581.85	7942.52	9033.44	9494.61	7469.87
7	61	8263.02	5756.22	7709.59	7991.12	8599.23	10386.9	9703.77	8570.38	8184.02	9533.22	9836.17	7973.57	5677.03	8479.69	8577.07	7940.11	8936.6	9400.46	7284.23

Plot ABC posterior marginal distributions

histogram(marginals[:ℓerr],title="ℓerr",labels="",size=(500,250),normalize=:pdf)

histogram(marginals[:ρ],title="ρ",labels="",size=(500,250),normalize=:pdf)

histogram(marginals[:αeff],title="αeff",labels="",size=(500,250),normalize=:pdf)

histogram(marginals[:βα],title="βα",labels="",size=(500,250),normalize=:pdf)

histogram(marginals[:βθ],title="βθ",labels="",size=(500,250),normalize=:pdf)

histogram(marginals[:γα],title="γα",labels="",size=(500,250),normalize=:pdf)

histogram(marginals[:γθ],title="γθ",labels="",size=(500,250),normalize=:pdf)

histogram2d(marginals[:γθ],marginals[:γα],labels="",xlabel="γθ",ylabel="γα",size=(500,250))

βμ = [mean(Weibull(marginals[:βα][i],marginals[:βθ][i])) for i=1:length(marginals[:βα])];
γμ = [mean(Weibull(marginals[:γα][i],marginals[:γθ][i])) for i=1:length(marginals[:γα])];

βσ = [std(Weibull(marginals[:βα][i],marginals[:βθ][i])) for i=1:length(marginals[:βα])];
γσ = [std(Weibull(marginals[:γα][i],marginals[:γθ][i])) for i=1:length(marginals[:γα])];

histogram2d(βμ,γμ,xlabel="Length of infectious period mean",ylabel="Length of recovery period mean",size=(500,325))

plt = palette(:default);
pltβγμ = histogram(βμ,normalize=:pdf,labels="",linecolor=:white,alpha=0.5,linewidth=0);
histogram!(γμ,normalize=:pdf,labels="",linecolor=:white,alpha=0.5,linewidth=0)
density!(βμ,c=plt[1],linewidth=3,linestyle=:solid,labels="contact interval β");
density!(γμ,c=plt[2],linewidth=3,linestyle=:dash,labels="infectious period γ");
plot!(xlabel="mean",ylabel="pdf",xlims=(12.25,15.0),
      size=(figwdth/2*λplot,fighght*λplot),
      tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
      titlefontsize=round(figtitlefontsize*λfont))

savefig("postmean.pdf")

# Master figure for all posteriors
# Ex of how to use blank plots to center and vary number of subplots in a plot
#pltblank = plot(legend=false,grid=false,foreground_color_subplot=:white);
#lay = @layout [_{0.165w} a{0.335w} b{0.335w} _{0.165w};c d e];
#plot(pltblank,plt1,plt2,pltblank,plt3,plt4,plt5,layout = lay,size=(figwdth*λplot,2*fighght*λplot))

plt1 = histogram(marginals[:ρ],title="ρ",labels="",size=(500,250),normalize=:pdf,ylabel="pdf",xlabel="",
                 xticks=0.0:0.025:0.05,linecolor=:white,linewidth=0)
plt2 = histogram(βμ,title="βμ",labels="",size=(500,250),normalize=:pdf,xlabel="",
                 xticks=12.4:0.3:13.0,linecolor=:white,linewidth=0)
plt3 = histogram(βσ,title="βσ",labels="",size=(500,250),normalize=:pdf,xlabel="",linecolor=:white,linewidth=0)
plt4 = histogram(marginals[:αeff],title="αeff",labels="",size=(500,250),normalize=:pdf,ylabel="pdf",xlabel="value",
                 linecolor=:white,linewidth=0)
plt5 = histogram(γμ,title="γμ",labels="",size=(500,250),normalize=:pdf,xlabel="value",
                 xticks=12.5:1.0:15.0,linecolor=:white,linewidth=0)
plt6 = histogram(γσ,title="γσ",labels="",size=(500,250),normalize=:pdf,xlabel="value",
                 xticks=1.8:0.4:3.0,linecolor=:white,linewidth=0)

plot!(plt1,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
           titlefontsize=round(figtitlefontsize*λfont))
plot!(plt2,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
           titlefontsize=round(figtitlefontsize*λfont))
plot!(plt3,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
           titlefontsize=round(figtitlefontsize*λfont))
plot!(plt4,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
           titlefontsize=round(figtitlefontsize*λfont))
plot!(plt5,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
           titlefontsize=round(figtitlefontsize*λfont))
plot!(plt6,tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
           titlefontsize=round(figtitlefontsize*λfont))
lay = @layout [a b c;d e f]
plot(plt1,plt2,plt3,plt4,plt5,plt6,layout = lay,size=(figwdth*λplot,2*fighght*λplot))

savefig("postprm.pdf");

$R_0$ posterior distributions

The $R_0$ value is given by $R_0 = \int^\infty_0 S_\gamma(t)\beta(t)dt$, where $S_\gamma(t)$ is the survival function of $\gamma$ and $\beta$ is the contact interval hazard. For a Weibull distribution, the survival function is $S_\gamma(t) = e^{-(t/\gamma_\theta)^{\gamma_\alpha}}$. Combining with the definition of $\beta$ we obtain
$R_0 = \int^\infty_0 e^{-(t/\gamma_\theta)^{\gamma_\alpha}}\frac{\beta_\alpha}{\beta_\theta}\left(\frac{t}{\beta_\theta}\right)^{\beta_\alpha-1}dt.$
We now estimate this integral by truncating it to a finite domain with error at most tol, and then use a trapezoidal approximation with npts many sample points.

function R₀(β_α::Real,β_θ::Real,γ_α::Real,γ_θ::Real;tol::Real=1e-6,npts::Integer=20001)
   # Find L₀ where tail of integral is at most tol when the tail is bounded by exp(-s/2)
   # You want I<tol*β_θ^β_α*γ_α/β_α/γ_θ^β_α
   # d/dx F = exp(-s/2) for F = -2 exp(-s/2); so I = 2 exp(-s/2)
    L₀ = -2*log(0.5*tol*β_θ^β_α*γ_α/β_α/γ_θ^β_α);
    L₀ = L₀>1 ? L₀ : 1.0;
    
    # Find where in [⋅>2r], the exponential is above the poly and may dominate with exp(-s/2)
    r = β_α/γ_α-1;
    L₀ = L₀ > 2*r ? L₀ : 2*r;
    while r*log(L₀) > L₀/2
        L₀ *= 2;
    end
    
    # Compute definite integral on ∫^L₀_0 by trapezoidal rule
    #  Note since u^(β_α/γ_α-1) can blow up at 0 in an integrable fashion
    #  we truncate the origin sample point so an inessential singularity
    #  does not cause problems
    xs = LinRange(L₀/npts,L₀,npts-1); Δ = xs[2]-xs[1];
    ys = β_α*γ_θ^β_α/β_θ^β_α/γ_α*exp.(-xs).*(xs.^r);
    
    val = [Δ*(ys[i]+ys[i-1])/2 for i=2:length(xs)];
        
    return sum(val)
end;

marginals[:R₀] = [R₀(marginals[:βα][i],marginals[:βθ][i],
                     marginals[:γα][i],marginals[:γθ][i]) for i=1:length(marginals[:βα])];
pltR₀ = histogram(marginals[:R₀],normalize=:pdf,labels="",size=(450,250),bins=15,linecolor=:white,linewidth=0,
                  alpha=0.5)
flag = marginals[:R₀].==Inf; marginals[:R₀]=marginals[:R₀][.!flag];
println("The mean R₀: ",sum(marginals[:R₀])/length(marginals[:R₀]))
println("A 95% CI for R₀: ",[quantile(marginals[:R₀],0.025),quantile(marginals[:R₀],0.975)])
density!(marginals[:R₀],c=plt[1],linewidth=3,labels="")
plot!(xlabel="R₀ value",ylabel="pdf",xticks=0.75:0.5:2.25,
      size=(figwdth/2*λplot,fighght*λplot),
      tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
      titlefontsize=round(figtitlefontsize*λfont))

The mean R₀: 1.2803727577313408
A 95% CI for R₀: [0.8880207220940958, 1.944244680190886]

savefig("postR0.pdf")

lay = @layout [a b];
plot!(pltR₀,ylabel="");
plot(pltβγμ,pltR₀,layout=lay,
      size=(figwdth*λplot,fighght*λplot),
      tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
      titlefontsize=round(figtitlefontsize*λfont))

savefig("mst_meanR0.pdf");

Parameter Posterior CI's

dftemp = DataFrame(:βα=>quantile(marginals[:βα],[0.0275,0.975]),
               :βθ=>quantile(marginals[:βθ],[0.0275,0.975]),
               :γα=>quantile(marginals[:γα],[0.0275,0.975]),
               :γθ=>quantile(marginals[:γθ],[0.0275,0.975]),
               :ρ=>quantile(marginals[:ρ],[0.0275,0.975]),
               :αeff=>quantile(marginals[:αeff],[0.0275,0.975]),
               :βμ=>quantile(βμ,[0.0275,0.975]),
               :βσ=>quantile(βσ,[0.0275,0.975]),
               :γμ=>quantile(γμ,[0.0275,0.975]),
               :γσ=>quantile(γσ,[0.0275,0.975])
                )

2 rows × 10 columns

	βα	βθ	γα	γθ	ρ	αeff	βμ	βσ	γμ	γσ
	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	5.05714	13.5036	4.97674	13.4295	0.00792989	0.706579	12.4068	1.61697	12.5304	1.93351
2	9.54716	13.5076	7.9382	15.0333	0.0495332	0.992436	12.8249	2.81312	14.1182	2.91921

Pearson correlation coefficients

μs = Dict{Symbol,Float64}(); σs = Dict{Symbol,Float64}();
for key in keys(marginals)
    μs[key] = sum(marginals[key])/length(marginals[key]);
    σs[key] = √( sum((marginals[key].-μs[key]).^2)/length(marginals[key]) );
end
n=length(keys(marginals));
PC = Matrix{Float64}(undef,n,n); id₁=0;
for key₁ in keys(marginals)
    id₁+=1;
    id₂=0;
    for key₂ in keys(marginals)
        id₂+=1;
        PC[id₁,id₂] = sum( (marginals[key₁].-μs[key₁]).*(marginals[key₂].-μs[key₂]) )/(σs[key₁]*σs[key₂]);
        PC[id₁,id₂] *= 1/length(marginals[key₁]);
    end
end;

dfpc = DataFrame("prm"=>[key for key in keys(marginals)]);
pos = 0;
for key in keys(marginals)
    pos += 1;
    dfpc[!,string(key)] = convert(Vector,PC[:,pos]);
end
println("Pearson correlation coefficients:")
dfpc

Pearson correlation coefficients:

8 rows × 9 columns

	prm	R₀	ρ	ℓerr	βθ	αeff	βα	γα	γθ
	Symbol	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	R₀	1.0	0.593828	0.709237	0.453027	0.0258399	0.453027	-0.106132	0.792204
2	ρ	0.593828	1.0	0.0189204	0.349604	0.0488572	0.349573	-0.069226	0.44916
3	ℓerr	0.709237	0.0189204	1.0	0.383182	-0.00831572	0.383154	-0.0803331	0.592233
4	βθ	0.453027	0.349604	0.383182	1.0	0.0756349	1.0	0.124744	-0.067327
5	αeff	0.0258399	0.0488572	-0.00831572	0.0756349	1.0	0.0756013	0.0805348	-0.00315127
6	βα	0.453027	0.349573	0.383154	1.0	0.0756013	1.0	0.12476	-0.0673385
7	γα	-0.106132	-0.069226	-0.0803331	0.124744	0.0805348	0.12476	1.0	0.0897513
8	γθ	0.792204	0.44916	0.592233	-0.067327	-0.00315127	-0.0673385	0.0897513	1.0

Analyze Best Fit

Simulate the best fit

df[!,1] = Symbol.(df[!,1])

# Find best fit
idℓerr = 1;
while df[idℓerr,1]!=:ℓerr
    idℓerr+=1;
end

ℓs = [df[idℓerr,k] for k=2:ncol(df)]
ℓbf = minimum(ℓs);
pos = findfirst(ℓs.==ℓbf) + 1;

# Prepare the data vector
vkeys = df[:,1]; prm,_=rdprm(df[:,pos],vkeys);
data!(prm);

println("The best fit parameters:")
for key in keys(prm)
    println(key,": ",prm[key])
end

The best fit parameters:
αeff: [0.9698334727742906]
nndsmp: [2500.0]
αL: [14.0]
nnd: [2500.0]
dwnsmp: [1.0]
βθ: [13.503792220531302]
atol: [1.0e-5]
Δtmin: [1.0e-5]
ℓerr: [746.6170971719679]
γα: [6.052065267527739]
γθ: [13.617999656279881]
fˢη: [3.301535733241945e-5]
yˢrg₀: [0.0, 36500.0]
yᵛrg₀: [0.0, 0.0]
λ: [1.5]
βα: [5.238644772185896]
yⁱrg: [0.0, 22.0]
ρ: [0.04476143286480499]
nΔtfail: [NaN]
rtol: [0.001]
yᵛrg: [0.0, 22.0]
yˢrg: [0.0, 36565.0]
fⁱη: [0.07407407402502883]
yⁱrg₀: [0.0, 14.0]
Trg: [0.0, 65.0]

# Run the simulation
ysol,yʳsol = pdesolve(;prm=prm);

Simulation progress: 0.015625/1
Simulation progress: 0.03125/1
Simulation progress: 0.046875/1
Simulation progress: 0.0625/1
Simulation progress: 0.078125/1
Simulation progress: 0.09375/1
Simulation progress: 0.109375/1
Simulation progress: 0.125/1
Simulation progress: 0.140625/1
Simulation progress: 0.15625/1
Simulation progress: 0.171875/1
Simulation progress: 0.1875/1
Simulation progress: 0.203125/1
Simulation progress: 0.21875/1
Simulation progress: 0.234375/1
Simulation progress: 0.25/1
Simulation progress: 0.265625/1
Simulation progress: 0.28125/1
Simulation progress: 0.296875/1
Simulation progress: 0.3125/1
Simulation progress: 0.328125/1
Simulation progress: 0.34375/1
Simulation progress: 0.359375/1
Simulation progress: 0.375/1
Simulation progress: 0.390625/1
Simulation progress: 0.40625/1
Simulation progress: 0.421875/1
Simulation progress: 0.4375/1
Simulation progress: 0.453125/1
Simulation progress: 0.46875/1
Simulation progress: 0.484375/1
Simulation progress: 0.5/1
Simulation progress: 0.515625/1
Simulation progress: 0.53125/1
Simulation progress: 0.546875/1
Simulation progress: 0.5625/1
Simulation progress: 0.578125/1
Simulation progress: 0.59375/1
Simulation progress: 0.609375/1
Simulation progress: 0.625/1
Simulation progress: 0.640625/1
Simulation progress: 0.65625/1
Simulation progress: 0.671875/1
Simulation progress: 0.6875/1
Simulation progress: 0.703125/1
Simulation progress: 0.71875/1
Simulation progress: 0.734375/1
Simulation progress: 0.75/1
Simulation progress: 0.765625/1
Simulation progress: 0.78125/1
Simulation progress: 0.796875/1
Simulation progress: 0.8125/1
Simulation progress: 0.828125/1
Simulation progress: 0.84375/1
Simulation progress: 0.859375/1
Simulation progress: 0.875/1
Simulation progress: 0.890625/1
Simulation progress: 0.90625/1
Simulation progress: 0.921875/1
Simulation progress: 0.9375/1
Simulation progress: 0.953125/1
Simulation progress: 0.96875/1
Simulation progress: 0.984375/1
Simulation progress: 1.0/1
Number of times Δt was too small for refinement: 0

Plot equation coefficients

plot(:α;prm=prm)

plot(:β;prm=prm)

plot(:γ;prm=prm)

plot(:Weibull;prm=prm)

plot(:fˢ;prm=prm)

plot(:fⁱ;prm=prm)

Plot best fit solution

plot(ysol)

plot(ysol,yʳsol;prm=prm)

plotbd(ysol;prm=prm)

Plot errors with ODH

df_yⁱ = CSV.read("ODH_snipdaily.csv",DataFrame);
first(df_yⁱ,3)

3 rows × 12 columns

	time	0-9	10-19	20-29	30-39	40-49	50-59	60-69	70-79	80+	daily_confirm	cum_confirm
	Date	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64
1	2020-11-15	559.202	578.048	1628.5	1341.0	1335.75	1421.5	1134.5	679.25	474.75	9152.5	3.60476e5
2	2020-11-16	561.538	580.462	1642.2	1354.4	1344.0	1432.4	1138.4	676.2	486.2	9215.8	365386.0
3	2020-11-17	556.457	575.21	1641.17	1368.33	1363.33	1440.5	1154.5	685.167	499.667	9284.33	370264.0

# Extract error for this trajectory
rms,yⁱ_daily = ℓerr(ysol;prm=prm);

# Compute 95% quantiles
qlow = [quantile(dftrjabc[i,2:end],0.0001) for i=1:nrow(dftrj)];
qhgh = [quantile(dftrjabc[i,2:end],0.9999) for i=1:nrow(dftrj)];

npts = nrow(df_yⁱ);
#plot(df_yⁱ[!,:time],df_yⁱ[!,:daily_confirm],linewidth=1,labels="ODH raw",linestyle=:dash,
#     title="ODH Daily Cases",xlabel="date",ylabel="count",size=(400,225));
odhavg = Vector{Float64}(undef,npts);
@inbounds for i=1:npts
    Σ=0.0; ct = 0;
    for k=-3:3
        try
            Σ += df_yⁱ[i+k,:daily_confirm];
            ct += 1;
        catch
            true;
        end
    end
    odhavg[i] = Σ/ct;
end

plt = palette(:default);
plot(df_yⁱ[!,:time],yⁱ_daily[1:npts],linewidth=0,labels="",c=plt[1],
                     ribbon=(yⁱ_daily[1:npts]-qlow,qhgh-yⁱ_daily[1:npts]),
                     fillalpha=0.25)
plot!(df_yⁱ[!,:time],yⁱ_daily[1:npts],linewidth=3,labels="model",c=plt[1],
                     #ribbon=(yⁱ_daily[1:npts]-qlow,qhgh-yⁱ_daily[1:npts]),
                     xlabel="date",ylabel="daily cases")
plot!(df_yⁱ[!,:time],odhavg,linewidth=3,labels="ODH 7 day mvavg",c=plt[2],
                     linestyle=:dash);
plot!(legend=:bottomright)
plot!(size=(figwdth*0.75*λplot,fighght*λplot),
      tickfontsize=round(figtickfontsize*λfont),guidefontsize=round(figguidefontsize*λfont),
      titlefontsize=round(figtitlefontsize*λfont))

savefig("abcvar.pdf");